package Integration;

import MathLib.triangulation.SaveStructures.Node;

/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

/**
 *
 * @author mark
 */
public class GaussQuadrature2D5Point {
    protected double x[];
    protected double w[];

    public GaussQuadrature2D5Point() {
        x = new double[5];
        w = new double[5];
        x[0] = 0.0;
        x[1] = -Math.sqrt(245 - 14 * Math.sqrt(70)) / 21.0;
        x[2] = Math.sqrt(245 - 14 * Math.sqrt(70)) / 21.0;
        x[3] = -Math.sqrt(245 + 14 * Math.sqrt(70)) / 21.0;
        x[4] = Math.sqrt(245 + 14 * Math.sqrt(70)) / 21.0;

        w[0] = 128.0 / 225.0;
        w[1] = (322.0 + 13.0 * Math.sqrt(70)) / 900.0;
        w[2] = (322.0 + 13.0 * Math.sqrt(70)) / 900.0;
        w[3] = (322.0 - 13.0 * Math.sqrt(70)) / 900.0;
        w[4] = (322.0 - 13.0 * Math.sqrt(70)) / 900.0;
    }

    protected double gau(double ax, double bx,double ay, double by, IFunction2D func) {
        double res = 0.0;
        for(int j=0;j<5;j++){
            double temp=0.0;
            for (int i = 0; i < 5; i++) {
                temp += w[i] * func.calculate((ax + bx) / 2 + (bx - ax) / 2 * x[i],(ay + by) / 2 + (by - ay) / 2 * x[j]);
            }
            res += temp*w[j];
        }
        res *= (bx - ax)*(by - ay) / 4.0;
        return res;
    }
    protected double ga(double ax, double bx,double ay, double by,int n, IFunction2D func) {
        double res = 0.0;
        double hx = (bx-ax)/(double)n;
        double hy = (by-ay)/(double)n;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                res+=gau(ax+i*hx, ax+(i+1)*hx,ay+j*hy, ay+(j+1)*hy, func);
            }
        }
        return res;
    }

    public double Gauss(double ax, double bx,double ay, double by, double eps, IFunction2D func){
        int n=1;
        double prev = ga(ax, bx, ay, by, n, func);
        double next = ga(ax, bx, ay, by, n * 2, func);
        while (Math.abs(prev - next) > eps){
            prev = next;
            n *= 2;
            next = ga(ax, bx, ay, by, n * 2, func);
        }
        return next;
    }

    private double GaussOnTriangle(Node A, Node B, Node C, double eps, IFunction2D func){
        int n=1;
        FuncWithJacobian ff = new FuncWithJacobian(A, B, C, func);
        double prev = ga2(0, 1, 0, 1, n, ff);
        double next = ga2(0, 1, 0, 1, n * 2, ff);
        while (Math.abs(prev - next) > eps){
            prev = next;
            n *= 2;
            next = ga(0, 1, 0, 1, n * 2, ff);
        }
        return next;
    }

    private double ga2(double ax, double bx,double ay, double by,int n, IFunction2D func) {
        double res = 0.0;
        double hx = (bx-ax)/(double)n;
        
        for (int i = 0; i < n; i++) {    
            double hy = (by-ay)/(double)n;
            for (int j = 0; j < n; j++) {
                res+=gau2(ax+i*hx, ax+(i+1)*hx,ay+j*hy, ay+(j+1)*hy, func);
            }
        }
        return res;
    }

    private double gau2(double ax, double bx,double ay, double by, IFunction2D func) {
        double res = 0.0;

        for(int i=0;i<5;i++){
            double temp=0.0;
            by=0.5 -  x[i]/2.0;
            for (int j = 0; j < 5; j++) {
                temp += w[j] * func.calculate((ax + bx) / 2 + (bx - ax) / 2 * x[i],  (ay + by) / 2 + (by - ay) / 2 * x[j]);
            }

            temp *= (by - ay) / 2.0;
            res += temp*w[i];
        }
        res *= (bx - ax) / 2.0;
        return res;
    }

     public double g(double ax, double bx,double ay, double by, IFunction2D func,Node A, Node B, Node C) {
        double res = 0.0;
        func = new FuncWithJacobian(A, B, C, func);


        for(int i=0;i<5;i++){
            double temp=0.0;
            by=0.5 -  x[i]/2.0;
            for (int j = 0; j < 5; j++) {
                temp += w[j] * func.calculate((ax + bx) / 2 + (bx - ax) / 2 * x[i],  (ay + by) / 2 + (by - ay) / 2 * x[j]);
            }

            temp *= (by - ay) / 2.0;
            res += temp*w[i];
        }
        res *= (bx - ax) / 2.0;
        return res;
    }
}
